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Abstract 

This document reports investigations of models of multiphase flows us- 
ing Lattice Boltzmann methods. The emphasis is on deriving by Chapman- 
Enskog techniques the corresponding macroscopic equations. The singular 
interface (Young-Laplace-Gauss) model is described briefly, with a discussion 
of its limitations. The diffuse interface theory is discussed in more detail, 
and shown to lead to the singular interface model in the proper asymptotic 
limit. The Lattice Boltzmann method is presented in its simplest form appro- 
priate for an ideal gas. Four different Lattice Boltzmann models for non-ideal 
(multi-phase) isothermal flows are then presented in detail, and the result- 
ing macroscopic equations derived. Partly in contradiction with the published 
literature, it is found that only one of the models gives physically fully ac- 
ceptable equations. The form of the equation of state for a multiphase system 
in the density interval above the coexistance line determines surface tension 
and interface thickness in the diffuse interface theory. The use of this rela- 
tion for optimizing a numerical model is discussed. The extension of Lattice 
Boltzmann methods to the non-isothermal situation is discussed. 
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Preface 



This document reports on a project carried out by the authors in March-August, 
2000. While waiting for renewed funding the results where not properly written up, 
but left in the form found here. It now being clear that the project has definitively 
been discontinued, the unfinished report is deposited as e-Print as is, in the hope 
that some of it may be of some use to someone somewhere. 

As a short guide, the paper contains first a survey of the diffuse interface models, 
also known as Cahn-Hilliard models. None of this is new, but is here collected in 
one place. It is pointed out that two different versions of the diffuse interface theory 
are used in the literature. The paper thereafter contains an introduction the Lattice 
Boltzmann equations along standard lines, and then a rederivation of the model of 
Swift et al. Again, nothing here is new, but the calculations are spelled out in more 
detail than in Swift et al. 

The paper then contains a survey of other proposed LBE models (isothermal), chap- 
ters 5-7, where we reach different results than in the published literature. In fact, 
we find that one such model not only does not give the correct momentum equa- 
tion, but even displays a spurious density diffusion term in the continuity equation! 
On the other hand, we find that one other model presented in the literature actually 
gives the correct continuity and momentum equations. 

Section 8 contains a discussion of thermal multiphase Lattive Boltzmann equations. 
The investigations reported are not conclusive, and probably along the wrong path. 
The general setting, and the references, could however be of use. We believe that 
constructing thermal multiphase Lattice Boltzmann equations is a difficult problem. 
We doubt it can be done in the standard simplifying single relaxation-time approxi- 
mation. 

Section 9 contains a discussion of optimization issues for the model described in 
section 7. By changing the equation of state in the interface region, one may change 
any of either three properties, keeping the others constant. The properties are the 
characteristics of the stable phases, the surface tension and the interface thickness. 
If, for instance, one wants to simulate an interface on a given grid, then the interface 
thickness in the model must be significantly larger. Hence, by changing only inter- 
face thickness, one could simulate multiphase fluids with given density contrast and 
surface tension on different grids. The numerical results do not quite illustrate this 
fact, which we nevertheless believe should be important in practical use. 
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The work in sections 4 and 5 were done in collaboration with Massimo Vergassola, 
Observatoire de Nice, to whom we direct our heartful thanks for patient instruction 
and very useful advice on the calculations. 
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1 Introduction 



Multiphase flows are difficult because they involve thermodynamics (co-existing 
phases) kinetics (nucleation, phase transitions), and hydrodynamics (inertial ef- 
fects). The nature of a liquid- vapour interface is still partly an open problem from 
the physical point of view. Mathematical descriptions of a multiphase flow have 
hence the status of models, and can be classified as i) singular interface models, 
which go back to Young, Laplace and Gauss; ii) diffuse interface models, which 
were first developed by Maxwell, Gibbs, Van der Waals [|4^] and Korteweg [p3l], 
and in a more modern setting by Cahn & Hilliard [^, and iii) fully detailed statisti- 
cal mechanics models, which are still under development [O]. 



In numerical work the most common are methods based on the first kind of models, 
mainly because they do not require the resolution of the interface, and can therefore 
be simulated with less effort. We begin below with a brief discussion of such mod- 
els, and in which situations they are likely to be inadequate. We then turn to the 
diffuse interface theory, from which one can derive the singular interface theory as 
an asymptotic limit. 

If one is only interested in phenomena where the singular interface model is ex- 
pected to be correct, the diffuse interface theory can, from the computational point 
of view, be considered a method of solving a moving boundary problem by solv- 
ing a corresponding set of reaction-diffusion equations. The difference between a 
purely numerical procedure and the diffuse interface theory is that the latter is also 
a continuum thermodynamic model of the interface, and can therefore be compared 
with a wider set of data. 



Lattice Boltzmann methods are computational schemes to solve macroscopic equa- 
tions by solving the equations of a microscopic world, which is described by the 
desired macroscopic equations in an asymptotic limit. An review of early work in 
the area is [^, for a comparatively recent review, see [[T0|]. If the microscopic world 
would be the real world, we would have a method akin to Molecular Dynamics, 
inevitably an expensive and cumbersome way to solve the macroscopic equations. 
What makes Lattice Boltzmann an interesting method is that we have freedom in 
choice of the microscopic model. In fact, we can make a rather drastic departure 
from the real world, and consider a microworld where velocities only take values in 
a discrete set. If such a model is correctly set up, the resulting macroscopic equa- 
tions will still be the desired ones. 
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From the PDE standpoint, Lattice Boltzmann methods would seem to be waste- 
ful, since they require resolution of scales smaller than hydrodynamic. In fact, 
however. Lattice Boltzmann methods do not require resolution of very small scales, 
and should rather be thought of as finite difference schemes with auxiliary vari- 
ables. In this regard Lattice Boltzmann methods are superior to Lattice Gas models, 
which use Boolean variables, and require an additional averaging, and hence a sig- 
nificantly larger scale separation. The fluctuations on the other hand make Lattice 
Gases advantageous for some thermodynamic stability problems, such as those re- 
cently investigated by Boghosian and coworkers [^. 

The competitiveness of Lattice Boltzmann methods with standard PDE methods - 
at present, and in the future - will mostly be determined by how large an overhead 
the extra variables in Lattice Boltzmann represent, and by the overall versatility and 
stability of the method. We will not attempt such a comparison in this report. A fair 
summary of the literature is that the opinion is still divided, but that Lattice Boltz- 
mann methods have so far mainly been used in the physics community, where the 
interpretation in terms of a fictious microworld is considered more of an advantage 
than a draw-back. A related point of view is that Lattice Boltzmann methods are 
closer to the physics, in the same sense as kinetics is compared with hydrodynamics. 
In flows with complicated physics the hydrodynamical descriptions may not have 
been worked out, or may be unfamiliar, and in both cases using Lattice Boltzmann 
methods could then be a more effective solution. The example usually given in this 
context is flow in porous media [|T^ . 

Given the reasons why Lattice Boltzmann methods for multiphase could be inter- 
esting, it still remains to construct such models. In this report we will go through 
models proposed in the literature, and find them lacking in several respects. Some 
of this work, presented below in sections 4-6 repeats calculations in the literature, 
while some of it is new. In section 7 we show that one of the schemes proposed in 
the literature on the other hand gives correct continuity and momentum equations 
for an isothermal multiphase flow. The construction involves an auxiliary quantity, 
and an expansion which mixes orders for this auxiliary quantity. In section 8 we 
show that the model of section 7 can be generalized to give non-isothermal multi- 
phase flows. We note that other Lattice Boltzmann models for non-isothermal flows 
have been introduced in the literature, but do not review them in detail. In section 9 
we discuss optimizing issues by modifying the free energy density in the interface 
region. 



5 



2 Models of multiphase flows 



2.1 The singular interface models 

In the singular interface model a surface is prescribed separating the two phases. 
The free energy includes a term proportional to the area of the surface, and spatially 
located with support on the surface, hence the name singular interface model. 

For definitiveness, assume that the fluid obeys Navier-Stokes equations in both 
phases, and that the interface is an impermeable separating surface. A necessary 
boundary condition is then that the normal velocity components are continuous 
across the interface 

V ■ = V ■ n\+ = Vn (2.1) 
If the fluids are viscous we also have that the tangential components are continuous: 

{v — {v ■ n)n) \_ = {v — {v ■ n)n) |_|_ (2.2) 

If the velocity gradients at the interface vanish, so that the interface and its surround- 
ings are locally at rest, and if the surface tension does not depend on the position 
along the surface, then there is a pressure jump across the interface according to 
Laplace's law 

P+ - P_ = 2aK, (2.3) 

where a is the surface tension and K, is the mean curvature of the surface, oriented 
such that the + side is on the inside and the — side is on the outside. If the fluid is 
not at rest, to the right hand side of ( |2.3| ) should be added the difference of viscous 
stresses, and if a depends on the position on the interface, e.g. through a depen- 
dence on temperature, there is a further term involving the gradient of a (see [^, 
§ 60). In general, singular interface models lead to free boundary problems with a 
stress balance at the boundary, and can include both heat and mass flow through the 



interface [14 1. 



Singular interface models are likely to be incorrect if applied to processes on length 
scales similar to the thickness of the interface. The standard list of such phenomena 
is (see e.g. [@]): 

• motion near a critical point, where the interface thickness is large; 

• motion of a contact line between a liquid and its vapour along a solid wall, 
see also [^; 
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• changes in topology, such as coalescence of liquid droplets or breakup of 
bubbles; 

• spontaneous formation of bubbles in overheated liquids (nucleation); 

• spontaneous separation into thermodynamically stable phases after a rapid 
quench from a higher temperature, where only one phase is stable (spinodal 
decomposition). 

The first three of these examples are from near-equilibrium statistical mechanics, 
while in the last two inertial effects are weak, at least initially. The physics of sys- 
tems that are both far from equilibrium and coupled to macroscopic transport is 
not systematically known, but is likely to give rise to new phenomena. One exam- 
ple could be bubble formation (nucleation) in the presence of strong temperature 
gradients, as in a liquid close to a heated wall. 



2.2 The diffuse interface theory: statics 

In modem terms, the diffuse interface theory is a Landau mean-field theory. In the 
isothermal setting, as considered by Van der Waals and Cahn-Hilliard, the (Helmholtz) 
free energy functional is assumed to depend on density and density gradients as 

J^= j^pfip,T) + ^K\Vp\'dV (2.4) 

where f{p,T) is the bulk free energy density (per unit mass) of a homogeneous 
phase with density p at temperature T, and K is a function of p and T. For simplic- 
ity we will in the following discussion take K a constant. 

A thermodynamically stable state at fixed mass is a minimum of the grand potential 
A = J-' — pM, where p is the chemical potential per unit mass. If there is a single 
spatially homogeneous global minimum at p, the free energy is 



T = Vpf{p,T) (homogeneous phase) (2.5) 

The thermodynamic pressure is P = — |^|(m,t) (derivative at constant mass and 
temperature). The chemical potential per unit mass is related to the free energy by 
p = ^\{v,T) (derivative at constant volume and temperature). Since M = pV, 
( P3D leads to 

P = p'f{p) ^i = f{p) + -pf\-p) (2.6) 
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Alternatively, if iJj = pf is the bulk free energy per unit volume, then 



P = p^'ip) - ^(p) p = ^'{p) (2.7) 

If, however, the grand potential has more than one global minimum, i.e. if there are 
two or more coexisting phases, then the mimimum of is found by a Euler-Lagrange 
variational equation: 

p = {pf)'p - KV^p (2.8) 
In one dimension, this equation can be integrated once to give 

{pf-p^i) = \K\d,p?-P (2.9) 

We expect there to be solutions of these equations which consist of domains where 
p is almost constant, and equal to one of the minima, separated by flat interfaces. 
Wherever V^p can be neglected, the Lagrange multiplier p is equal to pf'+f, which 
by ( [2.6[ ) is the chemical potential. If density varies but in one direction, then ( |2.9D is 
generally valid, and when then | Vpp can be neglected, the constant of integration P 
on the right hand side of is by ( ^!6| ) equal to the pressure. The thermodynamic rules 
that the chemical potential and the pressure in two coexisting phases in equilibrium 
are equal, can then be seen as conditions that the solutions of ( |2.8| ) and ( |2.9| ) are of 
the desired type. 

If we multiply with the gradient of p, we see that it is equivalent to 

d^T^f, = T^^ = p{f-p)S^^-K(d^pdpp-^\Vp\'5^^'^ (2.10) 

where T is a capillary stress tensor. In one dimension T has only one component, 
and (Og) reduces to (Ob. 



2.3 The diffuse interface theory: dynamics 

A difuse interface model of a multiphase flow is a set of equations for density (p), 
momentum density (pua) and energy density per unit mass (e), that reduces to the 
static case just discussed, if velocity is zero. These equations will be of the form 

dtp + dc,{pu^) = (2.11) 
pDtUf, = 9,(T,^ + ai'2) (2.12) 

pDte = da{{Tap + (T^^hup) - daqa (2.13) 
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where Tap is the reversible part of a stress tensor, a^'^j the hydrodynamic viscous 
stress tensor, and q an energy flux. Using the first (continuity) equation the left hand 
side of the two others can be rewritten dt{pUa) + dp{pUaUi3) and dt{p£)+da{p£Ua). 
The energy density consists of two parts: a kinetic energy and an internal 

energy density pei. The equations ( |2.1 1| ), ( |2.12| ) ( |2.13| ) must be supplemented with 



an equation of state, and a constitutive relation for the energy flux, which according 
to [0] contains a heat flux, and a "non-classical" term of interstitial working: 

Qa = -KdaT + Kp{V-u)daP (2.14) 

The appearence of the second term in the energy flux (but not in the entropy flux) 
can be traced back to an assumption by Van der Waals, that the gradient term in 
the free energy in fact stems entirely from the internal energy. It allows for some 
simplifications, in that the equations for the internal energy density ej and entropy 
density s read 

pDtei = da{t^daT) - p{V ■ u) + a^^ldaup (2.15) 
pTDtS = daiKdaT) + a^^^daUf, (2.16) 

where p is a pressure, related to density and temperature by the equation of state 
(p3|). If f[p, T, Vp] is the free energy functional density in ( |Z4| ) of a system at rest, 
ei[p, T, Vp] the internal energy functional density, and s(p, T) the entropy density , 



then f = Ei — Ts. The equation of heat transport (equation (2.16)) is therefore not 



independent, but can be derived from ( |2.15| ), the equation of state and thermody- 



namic identities, see § 49. The recent review [g] lists a number of problems in 
which diffuse interface models have been used successfully, for other recent papers, 
see [0] and [|21]]. 



Different dynamic diffuse interface models differ by the choice of q and the stress 



tensor T. Gibbs et al [48], who credit Lovett and Lebowitz & Percus, use for T the 



stress tensor in ( |2.1U| ). This choice has the perhaps undesired feature that it is not 



entirely local, since it depends on the chemical potential p: it assumes that some 
part of the fluid is at thermal equilibrium and at rest, so that p can be measured 
there. Another choice is to use equation ( |Z!8| ) to solve for p, and substitute that 
expression in ( |2.1(J| ): this is the form used in e.g. [Q]. We note that equation ( [2.8[ ) 



holds for a fluid at rest, but not necessarily for a fluid in motion, compare ([2.12[). 



We will not discuss further which of the most common forms in the literature gives 
the correct description of the dynamics of a diffuse interface. As a historic aside, 
it is however interesting to note that Korteweg in his 1901 paper P3|] derives by 
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symmetry arguments the following more general form of the inviscid stress tensor 
([p3|], p. 12, equation (20) includes also the viscous terms) 



+ hVapV^p+ dVlpP (2.17) 

with four undetermined coefficient functions (a, h, c and d) to describe the capillary 
forces. Korteweg's stress tensor reduces to the one of if one takes p equal to 
p{f ~ lA^^ equal to —K, a equal to K/2, and c = = 0. Korteweg further derives 
a surface tension (compare ( [2.22D 

which shows that the diagonal components of the stress tensor do not show up in 
the surface tension. 



We end this brief introduction by discussing the limits of validity of the Cahn- 
Hilliard free energy (p^. If is a constant, it follows that the surface tension 
decreases to zero along the critical line as cr ~ {Tc - TY with /X = |. The correct 
value of this critical exponent is however p f» 1.28 []5T|], which reflects the fact 
that a Landau theory does not correctly describe this second order phase transition. 
To match the observed dependence of the thickness one may take K a parametric 
function of temperature and density. It is still however an open question if such a 
modified Cahn-Hilliard theory correctly describes the interface close to a critical 
point. 



Along the coexistence line away from the critical point, the Cahn-Hilliard theory 
predicts a sharp interface, on the order of a few atomic diameters. This is in agree- 
ment with experiments, but nevertheless a conceptual problem, because if the inter- 
face is so sharp there is no reason why higher order terms should be relatively small. 
Since the singular interface theory can be derived from the diffuse interface theory, 
and since both have been applied to widely different flows, there should however be 
some meaning to truncating ( |2.4| ) after the first gradient term, perhaps as some sort 
of normal form. For a detailed discussion of gradient expansion resummation tech- 
niques of first principle statistical mechanics models, and of possible problems with 
the diffuse interface models which go outside the range of phenomena of interest to 
us here, see [ |I5| ] . 
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2.4 Derivation of the singular interface theory 

For completeness we derive here in the simplest case a singular interface model from a dif- 
fuse interface theory. 

If K is small and density gradients are close to zero everywhere but in a small part of 
the volume, then that part can be considered a fuzzy interface. Let us suppose that this in- 
terface passes through point P and is there approximated by a surface given parametrically 
hy z = ^(rix^ + r2y^). The two principal components of the curvature tensor (ri and 
r2) need not have the same sign. Let us further assume that density is only a function of 
z — i(rix^ + r2y^). A unit vector normal to the interface is formed by the density gradient 

n=^ = Jp£ff£$= (2.19) 

and points for small x and y approximately in the z direction. The stess tensor reads 

aai3 = (^-P + ^Pi^ Sai3 - Pin„n/3 (2.20) 

where P = p{fi — /) is the thermodynamic pressure in a homogeneous phase at density p, 

and pi is A'(|^)2, evaluated at a = z — ^(rix^ + r2y^). 

On the z axis the only non-zero component of the force field daC^ap is in the z direction: 

5QO"a/3(0, 0, z) = z - d^pi + pi(ri + r2)) (2.21) 

On the other hand, we know that at rest daCa/s = 0. If zq and zq are two points on the z 
axis on opposite sides of the interface we therefore have 

0= / do,aa(i{0,0,z)dxp = -{P{zQ)-P{zo)) + {n+r2) K\d,p\^dz (2.22) 

The pressure difference is proportional to the mean curvature C^^^^ ) and to the gradient 
contribution to the free energy in (2^). We can therefore identify that term as a free energy 
proportional to the area of the surface, i.e. a surface tension. 
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3 Lattice Boltzmann equation models in brief 



A physical system can be described on different levels, each of which correspond- 
ing to different resolutions in time and space. In statistical mechanics of essentially 
classical systems, such as ordinary liquids and gases, there are three important lev- 
els of description. The first (most detailed) level considers the motion of individual 
atoms and molecules. In a real liquid at room temperature and atmospheric pres- 
sure this implies spatial scale of about 10~^ m, and temporal scale of about 10~^^ 
s. Considerations at this scale are important e.g. in modelling chemical reactions in 
solution. 

The second level of description becomes valid when the actual motion can be sub- 
stituted with the mean number of particles to have momentum in a range ^ + A,^] 
at positions [x, x + Ax]. The condition for this level to be accurate is that fluctu- 
ations of this number are relatively small, or that the mean number is much larger 
than one. The system is then described by a chain of distribution functions, where 
X, t) is the probability density of finding a particle at position x with momen- 
tum ^ at time t, /2(^i, ^2- -^i- •^2- ^) is the joint probability of finding one particle at 
xi with momentum and another at X2 with momentum ^2, and so on. A descrip- 
tion of a gas or a liquid in terms of distribution functions is called kinetic theory. 
The simplest such description, valid for a very rare gas, is Boltzmann's equation, 
which in the single relaxation time approximation reads: 

^ + e-V/ = -^(/-r) (3.1) 

where f^'^ would be the distribution function of a liquid or a gas in termal equilib- 
rium with the same density, velocity and temperature. Boltzmann's equation can 
include a force field (acceleration term) and of course Boltzmann's full collision 
operator Vl, quadratic in /. In this paper we will however in accordance with most 
of the Lattice Boltzmann equation literature stay with the (simpler) single relaxation 
time approximation for fl. 

The first three moments of the distribution function / are respectively the density 
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(p), the momentum density, {fni) and the kinetic energy density per unit mass {sr)'- 

j f{lx,t)di= p{x,t) (3.2a) 
j U{1 f, t)di = p(f , t) u{x, t) (3.2b) 

/ ^•^^^'^'^^^^^ 

where m is the mass of a particle. We can introduce a local temperature field T 
such that Sk = \u'^ + ^RT, where D is dimension of space (2 or 3 in the cases of 
interest to us) and R the gas constant. The space integrals of density, momentum 
density and total energy density are integrals of the motion, namely total mass, total 
momentum and total energy. Large-scale modulations of these densities in space 
therefore have slow dynamics in time, compared to the other degrees of freedom 
on the kinetic level. A hydrodynamics of a physical system refers in general to a 
description only in terms of such slow, large-scale modes [||], which are therefore 
called hydrodynamic modes. Hydrodynamics in the ordinary sense, e.g. to describe 
the motion of water, is an example of such a description, albeit more complicated 
than for some other systems, since total energy is implicitly given by velocity, den- 
sity and temperature through the equation of state. 



A Lattice Boltzmann equation is a kinetic level description of a fictitious microworld 
with the same conserved quantities as the real world. If properly constructed (see 
below), it will then give rise to the same hydrodynamics. The smallest scale that 
needs to be resolved to compute the hydrodynamic equations in this manner is only 
such that hydrodynamics is valid for the scale of the object of interest. For prac- 
tical purposes. Lattice Boltzmann equations are therefore a sort of finite difference 
schemes, which contain auxiliary variables that stand for purely kinetic (not hydro- 
dynamic) modes. The number of auxiliary modes is typically 2-5 more than the 
number of hydrodynamic modes, which represents an overhead. 



A quick way to introduce the Lattice Boltzmann equations is to discretize space 
on a regular lattice, and momentum to a set of vectors in the dual lattice. Equation 
(PTTj) can then be solved to first-order accuracy in time 



Mx + e„t + At) - Mx,t) = --{f, - /f ) (3.3) 

T 

where the indices i label the set of momentum vectors, /, is the probability to be at 
lattice point x with momentum Cj at time t, and time At is such that a particle with 



13 



velocity Cj moves from x to x + Ci. The hydrodynamic variables are defined in a 
way similar to the continuum, i.e. 



^Wifi = p (3.4a) 

i 

^ WifiCi = pu (3.4b) 



\e,-u\^ = p^RT (3.4c) 



where the Wi are a set of weights. A given Lattice Boltzmann model is then de- 
termined by equation ( |J3D (or a generalization thereof), and how the equilibrium 
distribution function f^'^ depends on p, u and T. 

The a priori conditions for the same hydrodynamical equations to result from a 
Lattice Boltzmann model as from a real physical system, is that the symmetry of 
the lattice is sufficiently close to the full rotational symmetry of space. To capture 
correctly dissipative terms one needs that the invariant tensors form from the lattice 
vectors up to fourth order are isotropic: 

= 1 (3.5a) 

i 

WiCia = (3.5b) 

i 

^WiCiaCifS = T^^^5af3 (3.5c) 
i 

Wjeigeipej^ = (3.5d) 

i 

WieiaCipei^eie = T^'^\5ai3S-ye + Sa-yS^e + SaeSp-r) (3.5e) 

i 

Greek indices (a, 13...) here label spatial directions, i the lattice vectors and T*^^) and 
T*^''^ are lattice constants. Natural choices of lattices share with space the property 
that all odd order invariant tensors are zero, and the condition on second-order is 
easily satisfied. The technical difficulty overcome in the 80'ies (in the context of 
Lattice Gases) was to choose lattices with isotropic fourth order tensors, such as the 
hexagonal lattice in two dimensions. For faster convergence to the hydrodynamics 
it has sometimes been proposed to use lattices with isotropic sixth or higher order 
tensors [^0. 
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One presently popular model uses a lattice of 9 velocities in 2 dimensions, and 
is therefore known as the D2Q9 lattice. The velocities are the rest state (velocity 
eo = 0), along the vectors to the four nearest neighbours in a simple cubic lattice 
(velocities 61-64), and along the vectors to the four next-nearest neighbour (65-68). 
If lattice length is cAt, the velocities and weights are given in the following table: 



Direction 


Lattice vector 


Weight 


eo 


(0,0) 


Wo = 4/9 


ei 


6(1,0) 


wi = 1/9 


62 


e(0, 1) 


W2 = 1/9 


es 


e(-l,0) 


W3 = 1/9 


64 


6(0,-1) 


Wi = 1/9 


e5 


6(1,1) 


W5 = 1/36 


ee 


6(-l,l) 


= 1/36 


e? 


6(-l,-l) 


W7 = 1/36 


es 


6(1,-1) 


Ws = 1/36 



It is easily checked that all odd order invariant tensors formed from this lattice in- 



deed vanish, and that is satisfies (|3.5c[) and (p.5eD with T^^^ = ^ and T^^^^ 
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For the simplest applications it is sufficient to choose the equilibrium distribution 
function as 



1 + 



+ 



from which we can derive 

i 



T(2) 2TW 2T(2) 

^ p 



i 



T(2) 



(3.6) 

(3.7a) 
(3.7b) 
(3.7c) 

(3.7d) 



From ( [3.7cD it is natural to identify T*^^^p with a pressure term p. This means 
that the velocity of sound, Vg = y should be vT(^ (a constant). For the 
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D2Q9 models this constant is c/-\/3. The temperature RT is on the other hand 
^ifi'^i^ia — Ua){eia — Ua), a constant equal to T^^^. A model using 



and ( p.6p can therefore only describe isothermal hydrodynamics with an ideal gas 
equation of state, p = RTp. This is sufficient for deriving e.g. incompressible 
Navier-Stokes dynamics in a small Mach number expansion, but otherwise some- 
what limited. Two important issues in the following will be to generalize ( |33| ) or 



( p.6[ ) or both, to allow for flows which do not have the ideal gas equation of state. 



and which are not restricted to be isothermal. 

The usual presentation of Lattice Boltzmann methods (as above) uses a uniform 
grid in space. This limitation is not a principal one, by varying in an appropriate 
manner the weights Wi it is possible to use also non-uniform grids. These issues 
have not been investigated extensively in the literature, see however [ ]32| ] and [fT7|]. 
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4 Non-equilibrium thermodynamics model 



This model was developed by Swift and collaborators [M, B5p. The model has later 



been used in numerical work by several teams, for contributions from authors of 



this report, see [ pBQ and [pOfl. The basic idea is to assume that the distributions relax 
to a local equilibrium distribution function, which in a one-phase region coincides 
with the equilibrium in ordinary one-phase flows. In the interface region the local 
equilibrium is however assumed to have a pressure (stress) tensor as in the diffuse 



interface theory. The model uses the D2Q9 lattice, equation ( [3. 3D and assumes an 



(3.6) 



isothermal situation. The equilibrium distribution function is a generalization of 

The coefficients A, B, C, D, Ga,i3 and Aq and Co are chosen such that the zeroth 
and first weighted moments are p and pu, while the second moment is prescribed to 
be 

Wifl'^CiaCip = Pap + pUaUp (4.2) 

i 

where second order tensor P (in [M] referred to as a non-equilibrium pressure ten- 



sor) is chosen to be a stress tensor in a diffuse interface model: 

Pcpir) = pir)5ai3 + fidapd/sp (4.3) 
The diagonal component p(r) of the tensor P is 

pir) =po- KpV^p - ^|Vp|2 (4.4) 
and pq is the pressure in a Van der Waals equation of state 

Po = i—r - «P (4-5) 

1 — bp 

where T is temperature (assumed constant). The partial derivatives of the density 
are evaluated as discrete differences over neighbouring lattice points. In terms of a 
real system, this introduces a coupling between neighbouring sites, i.e. a potential 
energy. 

The kinetic energy density computed from the distribution function ( |4.1| ) would 
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in fact by ( ^T^ give rise to a local temperature field T(r) = y — i^V^p. The 
equilibrium distribution in this model is however computed not with this variable 
temperature, but with a global constant parameter. The use of the model is therefore 
restricted to phenomena where the back action of a temperature gradient on density 
and momentum can be neglected, basically a isothermal setting. 



The model starting from ( |4.3[ ) will give rise to Gallilean non-invarian terms in the 
hydrodynamic equations. It will become clear in the following that of the various 
discretization errors that appear in the Swift et al. model, there are three which are 
more serious, and should be considered of the same order as the diffusive/viscous 
terms. To the model will therefore be added a term Tap which eliminates two of 
these. 

Pafi{r) = p{r)5ap + i^dapd/sp + Tap (4.6) 



One result of of [45] is that the mass and momentum equations are: 



|^ + V.(pu)=0 (4.7) 

and 

^ + (u.V)(pu) = -Vpo + ^A(pu) + V.[A(p)V.(pu)] 

, (4.8) 
-5*^V.(V.(pu) + (pu).V) 

in which the kinematic viscosity v and bulk viscosity A of the fluid are determined 
by: 

_ 2r - 1 T(4) 



In fact, as we will see, equation ( |4.8D is only correct with the additional assumption 
that the fluid is close to incompressible, something which is also found in 

4.1 Method of successive approximations 

Following [Q] we will derive continuity equation and the Navier-Stokes equation 
by the method of successive approximations. The discussion in this section does 
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not contain other material than in p^, but is in places carried out in more detail. 
We start from equation ( |3.3| ) with a Taylor expansion of the l.h.s: 



k=l 



k=l 



with D = ((5t + Ciada)- We assume that the distribution function / is close to the 
equilibrium distribution function f^'* and expand in 5t- 

h = ft' + (/. - /D = fr-r {^{Df, + I^V.) + 0{6l) (4.10) 



Substituting ( |4-.10t ) to r.h.s of and retaining terms to order (5| 



2 



(4.11) 



So 



5. 



(5t + e.«a„)/r - 5* (r - (9i + eM'n' + (^(^i ) (4-12) 

By projecting on the hydrodynamic modes we will now derive the macroscopic 
equations to lowest and next lowest order in 5t- To obtain the continuity equation 
( p77| ) to lowest order, we sum over i and use ( [3.4aD and ( |3.4b| ), which gives 



= E = E ^^^'^ + E - ^(^*) 13a) 



Or 



5tP + daipUa) = 0{6t 



(4.13b) 
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To the next order in St'. 



(4.14) 



dt {dtp + dc,{pUc,)) + dc, dt{pUa) + 5/3 ^ WiGiaeipft'' 



(4.15) 



Substituting in eqn.( |4.15| ), eqn.( |4.13b| ) and ( |4.18| ) derived below we see that the 
0{5t) terms are actually 0{5l). The continuity equation is therefore: 



dtp + d^{pu^) = 0{5l 



(4.16) 



This means that there will not be any density diffusion terms in our equations, as 
indeed there should not be on general grounds . 



To obtain the momentum equation, ([4.8[), we multiply eqn.( [4.12] ) by Cj/j and sum 
over i: 



=dt ^ WiCipf-'^ + 5" X] '^i^ia^ipfr ' ^'(5, 
i i 

=dt{pUp) + da{PaP + pUaUfi) " 5f ( T - ^ 



(4.17) 
(4.18) 



+ 0{5l 



To lowest order we get the Euler equation which can be substituted back into the 
last parenthesis in eqn.( |4.15| ), to prove that the continuity equation holds to order 
5"^. Expanding P^p by equation ( P3| ) we have 



dt{pup) + daipUaUp) = -dppo + Kdp{pAp) + -9^|Vpp 



-da{ndapdpp) - daJ^aP + Ap + Bp + 



■/3 
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Equation ( |4.19| ) shows terms on the Euler level, where some terms have been simpli- 
fied using (|4.2|). We will now simplify further the three viscous-level terms denoted 
Ap, Bp and C^. Term Ap can be computed by substituting back eqn.( |4.19| ) 

dt{pup) = dt{dt{pup)) 

= -dt{d^{Pat3 + pUaUfs)) + 0{St) (4.20) 

In the following manipulations, we always imply all equalities to hold up to terms 
of order 5t . Therefore, 

Ai3 + =dt{da{PaP + pUaUp)) 

=dtd0o - dtd^ii^pAp) - dtdp (^IVpp) + 

dtda (Kdppdap) + dadtipUaUf^) + dadt^^p (4.21) 

We see that there are very many terms to order 5t- It will be useful to concentrate 
on the terms with smallest number of derivatives acting on the density. We will 
therefore, following Swift, neglect terms with as least as many derivatives as these: 

dtdfsinpAp) ^ 
dtd,[^\Vp\')^0 
dtda {ndppdap) ~ 

With this simplification, 

Af3 + Cf3 = dtdppo + da[Ui3dt{pUa)] + da[{pUc,)dtU,3] (4.22) 



ACl^ AC2p AC3^ 



These terms can be further simplified by 

ACl^ = dtdppo = dpj^idtp) = -dp (^d^{puS^ (4.23) 

AC2^ = -d^[uiid^{P^^ + puc^u^)] (4.24) 
Neglected furthermore these terms in equation ( 14.24D 

dady{KpAp) ^ 

dad^ (Kdgpdsp) ^ 

dadjT ^ 
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we have 

AC2p = -daiUpdaPo - UfsdjipUaU^)] (4.25) 

The term AC3 can be computed by using 

dt{pup) = updtp + pdtup = pdtup - Uad^{pu^) (4.26) 

So 

pdtUa = -dp{Pap + pu^up) + UadfsipUa) (4.27) 

Therefore 

ACSfs =dc,[ua{pdtUf3)] 

= - da{Uadf3Po - pU^U^d^Up] (4.28) 

Substituting equations ( |4.23| ), ( |4.25| ) and ( |4.28| ) back into ( |4.22| ) gives 

- da[updaPo + Maf?/3P0 " Upd^{pUaU^) - pu^u^d^up] (4.29) 



- 5/3 ( ^daipUa) ) - daiUfsdaPo + U^OpPo) + d^d^ipUa^UpU^) 



dp ( ^daipUa) ] - da 



dpo 
dp 



+ dad^{pUaUpU^) (4.30) 



The term Bp in the equation ( |4.18| ) can be computed by using the equilibrium dis- 



tribution function (4-. 1 ) 



dad^ ^^Y^Wif-'^eiaeipei^ =dad^ WiBugeigeipei^eie^ 

=dadjBT'^^\uj6ap + Up6a-y + UaSp^) (4.31) 



In eqn.( |4.31| ), the only free index is P, while a and 7 are summed over. The value 
of B given by Swift et al is p/T*^^)]] and therefore 



{2dpdaipUa) + ^{pUp)) 



(4.32) 



' The expression given by Swift et al is for the two-dimensional hexagonal lattice, for which 

T(2) = 3c2, = 3cV4 and hence BT^^^ = 
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Substituting now equations ( [4.30| ) and ( [4-.32D into ( |4.19| ) we have 



K 



6A2t - 1) 



T(2) 



pu 



VpVL 



r^^^-^da{pUa) - ^^daipUa) 



\{p)da{pUa) 



dp 



{Updap + Uadpp) + d^{pUaUpU^) 



(4.33) 



Some terms in the equation ( [4.33[ ) can be further simplified as: 



K 



nd/sipAp) + '-dfs (|Vpp) - Kdo,{d^pdpp) 

= Kdp{pAp) + Kdpp{dadpp) - K{dadpp)dpp - K{disp)Ap 
= Kpdp{Ap) 

So, the momentum equation is: 

dt{pUp) + daipUaUp) = -dpPo + uA^pUp) + df3[X{p)da{pUa)] 

dpo 



(4.34) 



-St T 



dp 



■da{Uadf3P + UpdaP) - daJ^afi + Kpdp{Ap) 



(4.35) 



The extra viscous terms in the momentum equation ( |4.35| ) are not Galilean invariant 
when density gradients are present. Comparing one sees that all of these terms can 



be derived from spatial derivatives acting on the pressure, i.e. term Cp in ( [4.18[ ). 
A possible interpretation is therefore that the Swift et al. condition on the sec- 
ond moments of the distribution function, i.e. equation (14. 6[), effectively involves 
an evaluation of a force term in Boltzmann's equation (compare section ^ below), 
which has to be done to better than linear order if unphysical viscous terms are to 
be avoided. 



Some of the non-Galilean invariant terms in ( }4.35D can be removed by an additional 
term JF^^ in the pressure tensor (equation (|4.6|)). First, we define: 



(4.36) 
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The equation ( [4-.35[ ) can then be rewritten: 



dt{pUp) + da{pUaUp) = - dppo + l^OaipdaU/s) + dp{\pdaUa) 

+ da[{jy- O^fsdaP - ^Uad(3p] (4.37) 
+ dp[XUadaP] + daTaP + Kpdp{^p) 

By a suitable choice of the term JF we can eliminate two of the three terms with one 
derivative of p 

^OLfi = ^{UfjdaP - Uadpp) - \u^d^p5af3 (4.38) 
and with this choice the momentum equation is becomes 

dt{pup) + daipUaUp) = - dppo + vda{pdaUp) + dp{\pd^Uc,) 

+ udaiufsdap) + Kpd^iAp) (4.39) 

With this modified model Swift et al. found qualitatively correct dynamics of a 
bubble in a liquid. This should be considered an experimental fact, difficult to 
explain theoretically. 

4.2 Summary 

The model of Swift et al. has some disadvantages: 

• The model is not Galilean invariant when the density gradients are present. 
By introducing the extra term JF, we can remove some of the non-Galilean 
invariant terms, but not all. 

• It is difficult to introduce the force terms into this model. We have to modify 
the continuity equation with a force terms F inside the derivative of pu. 

dtp + V.{pu + F) = (4.40) 

Therefore, this method can only be used if the force term F is negligibly 
small. 
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5 Luo model: LBE for non-ideal dense gases 



The model of Li-Shi Luo starts from Boltzmann equation for dense gases [^7]], also 
known as the Enskog equation: 



| + eV/ + a.V,/ 



n 



(5.1) 



where a and are the acceleration and collision operators. With the expansion of Q 
in a Taylor series, using a single relaxation time approximation, and assuming the 
fluid to be isothermal and incompressible, Luo obtains the following equation 



J=-rbpg{i-u).Wln{p'g) 



(5.2a) 
(5.2b) 



where h is the second virial coefficient of the equation of state for the hard-sphere 
system, and g is the radial distribution function. The equilibrium distribution func- 
tion is discussed by Luo in the continuum form 



eq 



P 



(27ri?T)^/2 
but is used in the discretized form ([3.6D. 



exp 



2RT 



(5.3) 



Solving (5.2) to first order in time we have 



{f{x,U)-r{x,^,t)) 



+ J{x,^,t)6t-a.V^f{x,^,t)6t (5.4) 



with a force term, a. V^/, written P7| ] 

a.V^/ = -puiO^r' [i^ - u) +eT'(eu)e] .a 



(5.5) 



The auxiliary variable is equal to a/RT (-^ in the D2Q9 model), and the weight 
weight uj(^) denotes the weight Wi for the momentum ^ = Cj. The force term J 
is proportional to the gradient of the function of density p'^g, which as in the Swift 
model is evaluated by discrete differences over lattice points. The interpretation is 
also that of a potential term. Equations ( |5.4D and ( |5.5[ ) together with (P^, ( |3.4aD , 



(3.4b), and a choice of a lattice, define Luo's Lattice Boltzmann model. 
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To proceed further we want to evaluate moments with respect to momentum of the 
two "forces" (the interior force J, and the exterior force through the acceleration a). 
The forcing term of eqn.(|5.5p is 



Fi = -Wip 



[Ci - u) + e,; 



T(2) 



.a 



and satisfies the following constraints: 



i 



T(2) 



^(wieiaei/3).a = -pa„ 



T(2) 

= piu^ap - upac,) 
The interior force term Jj satisfies the following constraints: 

i \ i i / 

= —bpg{pu — pu) .V ln{p^ g) = 



(5.6) 

(5.7a) 
(5.7b) 

.a 

(5.7c) 



(5.8a) 



^{JiCia) = -bpg I ^(/reiaCi/j) - ^{fl'^eio)u0 J dii\og{p^g) 

i \ i i / 

= -bT'^'a^p^a) (5.8b) 



-bpg ~ X^(/reiaei;3)M/3j .Vln{p^g) 

h[u^Uf,M.V - T(2)(m^c'^ + M^5«)](p2^) (5.8c) 



5.1 Chapman-Enskog analysis 

In this section we will derive the macroscopic equations with the Chapman-Enskog 
method. The earlier presented method of successive approximations can be consid- 
ered a special case of the more general and more flexible Chapman-Enskog method. 
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First, we apply the Taylor expansion to equation ( |5^ ) 

oo 

^ '-V^M^,t) = -^{M^,t) - /r(x,t)) + J,e - F,e (5.9) 

Tl. T 

ri=l 

where e = 5t and "D = ((9( + Ciada). We now introduce the following expansion: 

oo 



Ti=0 
oo 



dt=Y,^''dnt (5.10b) 
Tl=0 

a„ =9o„ (5.10c) 

The meaning of the expansion of / in powers of e is clear. The meaning of the for- 
mal expansion of the time derivatives in derivatives with respect to different "times" 
is that / is assumed to have support in the frequency domain on bands that are sep- 
arated by powers of e. We will look for an equation involving a "slow time", that 
is low-frequency components of /. Such components are Fourier transforms of /, 
where sufficienty fast oscillations are averaged out. A rigorous theory of these pro- 
cedures can be given in homogenization theory In general one can also assume 
there to be fast and slow spatial variables (small and large scales). For the present 
problem it is however enough to assume one spatial scale, hence ( |5.10cp 

With order accuracy, equation ( |5.9[ ) can be rewritten: 

i-(O) , ^.(1) , ^A-i) _ req ' I - t" _o \ „ T 



/. = fr + ^fr + ^'fr = /r - - + -v /, + -j,e - -F,e (5.11) 

^ V 2 y g g 
Group the right-hand side of the eqn.( |5.1 1] ) by the order of e. 

+ e-Ji-e-Fi (5.12) 
9 9 

To successive orders in e we find: 

/r (5.13a) 



e°: 















.(0) 
i 

'^9 9 9 



-Vtjr + -J,--F, (5.13b) 
9 9 9 
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The goal of the calculation is to project on the hydrodynamic modes of density 



and momentum. Equation ( |5.13aD shows that the zeroth order distribution function 
agrees with the equilbrium distribution function, which shares the same density and 
momentum. The higher order distribution functions therefore obey the following 
constraints: 



in) 







[n > 0) 



(5.14) 

$^(i^^"^e,«) = in>0) (5.15) 

i 

Summing equation ( |5.13b[ ) over i, using equations ( |5.7aD , and ( |5.8a| ), gives 



dtoP + doaipUa) = 



(5.16) 



To second order in e, we sum equation ( |5.13c[ ) over i, substituting the first order 
solution ( |5.13b| ) which gives 



+ (i-7)E»../^"-9'.E/^°' (5-17) 



or 



(5.18) 



Combining the first and the second order results for /j by dt = dt^^ + edt^ and 
recalling that e = 5*, we have a continuity equation with an (unphysical) density 
diffusion term: R 



dtp + da{puc 



St 



da {paa - bT^'^^daip'^g)) 



(5.19) 



The first order of the momentum equation can be produced by multiplying equation 
( |5.13b[ ) with Cifs, then summing over i: 



^The author of [ p7| ] claims there is no such density diffusion term, but we find that it is there. 
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By using equations ( |5.14| ), ( |5.8b| ), and ( |5.7b| ), we have 

dto{pUf3) + daipUaUfs) = -T(2)(9/3P - hT^'^'^dp{p^g) + pap 

= -dpP + pap 

where 

P = T(2)p(l + 6p^) 



(5.20) 



is a non-ideal pressure. Equation ( p.20D is hence Euler's equation. The second order 
momentum equation is 

i i i 

Subtituting the equations ( |5.15| ), ( ^J\) and ( |5^ ) into above we have 

9tApup) = - 1) ^-ni'j - ^^*o E^'^^^*/^ - ^^^^p) (5.21) 

where the first-order momentum flux tensor U^^^^ is f-^^eiaCip and therefore 



i/3 



= — do, E Ao iT^eiaCip +-da ^{JiCiaeip - Fidaeip) (5.22) 
where n^°j is the zeroth-order momentum flux tensor. Furthermore 

Q,/3 + pUaUp) + T^^^ [5q,(pMq,)5q,/3 + da{pUp) + (9/3(p-Ua)] 
= T(2) [dt^p + 9„(pMc,)](5«/3 + dt^XpUaUp) + T(2) [9c,(pM/3) + 9/3 (pM„)] 

=dt,{pUo,up) + T(2) [5,(pM^) + 9^(pM„)] + (5.23) 
Substituting back the above results into ( |5.21D , we have, 

dtiipUp) =T^^^ ^^^^ ^ da[da{pUp) + dp{pUa)] + ^^^^ ^ dadtoipUaUp) + 

^ dto "^{JiCip - Fidp) + i(9„ E( JieiaCj/j - FiCiaCip) (5.24) 
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Combining the first and the second order results, equations ( |5.20[ ) and ( |5.24D , for fi 
by dt = dto + edt^ and recalling that e = 5t, we then have the Navier-Stokes level 
momentum equation, 

dt{pup) + daipUaU/s) = -OpP + vda [da{pup) + dp{pUo)\ + pap + (5.25) 
where the kinematic viscosity is 



-5tin the D2QQ model 



2^? V % 

and various extra viscous-order terms are grouped together in £p: 

2r-g 



£(3 =St 



^9 



-dadtoipUaU/B) + 
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29 



-dp{p^9) - pap 



+ 



1 



^t-da {h[uc,up\i.y - T'^^\uadfi + Ufidc,)]{p9'^) + p{uaap + upa^)] (5.26) 

This momentum equation is not Galilean invariant when density gradients are present, 
and contains some unphysical terms. If we assume that Mach number is small (M 
formally of the same order as e in the expansion), then the extra term dadt^ {pUaUp) 
in 8p can be further simplified as: 

dadtoipUaUis) =da{Ua[dto{pUi3) - U/sdtoP] + M/3<9to (pUa) } 

=da{ — Uad^{pUj3U^) + UaUf3d-y{pU^) — U^d'j{pUaU.y) 

- {uadf^P + updaP) + p{uf3aa + u^ap) ] + 0{5t) (5.27) 

Terms of cubic order in Mach number {0{v?)) can now be neglected, and we sim- 
plify to 

dadtoipUaUp) = - T^'^^daiUadfSp + Upd^p) ~ VT^'^'^ dJ(lLa_dfi{p^ 9) + Ufido,{p^9)\ 

+ 5a [p(M/3a, + M„a^)] + 0{bt) + 0(u3) (5.28) 
which allows the extra terms in the momentum equation to be rewritten as with 

2t — 9 
•^9 
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Summary 

This model is correct to first order (inviscid equations), but incorrect to second 
order (viscous equations). The continuity equation contains density diffusion 
terms, and the momentum equation various unphysical terms, somewhat sim- 
ilar to the ones appearing in the model of Swift and collaborators. 
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6 Intermolecular interaction model 



This model was introduced by Qian, D'Humieres and Lallemand and further 
developped by He et al [|T9|]) and Chen et al. [[TTl]. The starting point is the Boltz- 
mann equation with a force field F in the single relaxation time approximation 
(compare (pl^) 



/-/ 
A 



eq 



(6.1) 



The force F is not specified. In the case at hand we want eventually to recover 
a diffuse interface theory, and F should therefore include density gradients. For 
different proposals, see [ [T9| ] and [|11]]. The ansatz made in is that the gradient 
Vg/ can be approximated by V^/'^'^. Since in the continuum we have a Maxwell- 
Boltzmann equilibrium distribution function. 



/ 



we further assume V^/ 
in the form 

dl 
dt 



eq 



(27r/?T)^/2 



exp 



2RT 



(6.2) 



'^wf^'^ and obtain an approximate Boltzmann equation 



f_feq F.(^-U) 



A 



RT 



f 



eq 



(6.3) 



By discretizing on a lattice and integrating over one time step to first order we obtain 
a Lattice Boltzmann scheme as 



fi{yi + eiSt,t + 6t) - /i(x,t) 



/(x,t)-r^(x,t) , F. (e,-u) 



+ 



(6.4) 



since in the isothermal models Lattice Boltzmann models RT = T^'^\ Comparing 
we see that Luo's model only differs from ( |6.4D by a definite choice of the force 
term F, and by a slightly different form of the external force term in (P3|). 

The continuity and Navier-Stokes equations will be derived by the method of suc- 
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cessive approximations. We note the constraints 



T(2) 

Ik. 

"Y(2) 



{pUa - UaP) = 



(6.5a) 



J2 ^^^^^ =^ ( ^io.^i?f? -^oY^ eipf-'^ j 

i \ i i / 



(6.5b) 



(6.5c) 

For the D2Q9 model the last expression simplifies to 



pFpu^ + pF^up -{F.\i)pu^ui3 



(6.6) 



We now sum the moments of the collision operator = — l/T(/j — f^'^) over i. 
The zeroth moment to first order reads 



dtP + da{pUa) = 0{St) 



(6.7) 



and the next 



i i i 



(6.8a) 



dt {dtp + daipUa)) + 9a dtipUa) + 5/3 ^i^^if^fr 



StrdaipFa) 



(6.8b) 
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Substituting in eqn.( |6.8bD , eqn.(|^) and ( |6.1CI| ) derived below, we see that, 

0{5t) = - 5i (^r - \0{5t) + 0{5t) + d^{pF^)] - Strdo^ipF^) 

with this, the continuity equation becomes. 



dtp + daipuc 



(6.9) 



To obtain the momentum equation, we multiply Vti by Cj/j and sum over i. To lowest 
order we get an Euler equation 



i i i 

dt{pUp) + <9„(pMaM/3) = -dfiPo + pFp + 0{6) 



(6.10) 



where Pq = T^'^^p is the ideal gas pressure. Pi given by pF^ = —d^Pi is a correc- 
tion, and the total pressure is P = Pq + Pi. We note that P can contain density 
gradient terms by a suitable choice of F. To next order in 5t we find 

dtipUfs) + da{pUaUfs) = - dfsPo + Fpp + 5t(^ - (<9t + ei^^do)'^ ^ f-'^en^ 



- hT{dlpFfi) + d^^pFpu^ + pF^up)\ + 0(5^) (6.11) 

In equation ( |6.1 1| ), a term (9q,(pMq,m^F.u) has been neglected because the order of 
this term is 0(u^), and F will eventually also have to be taken small. We have 



Y(4) 



=at[at(pu^) + dc,{pUaUp)] + 2T^'^^df3[dtp + <9c,(pM«)] + 



(6.12a) 
(6.12b) 



Y(4) 

dadtipUaUf}) + Y{i)dada{pUl3) + 2 



^(4) 



T(2) ^2 (pM,) (6.12c) 



For the D2Q9 model T*^'') = (T*^^))^, and the last parenthesis in above vanishes. By 
substituting the continuity equation ( p^ ) and momentum equation ( |6.1(JD , with the 
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accuracy to order of 5t, into ( |6.12cD we have 



A(3 = - T^'^'dtdpp + dt{pFp) + dadtipu^Ufs) + ^(^^a^alpM/s) + (6.12d) 



Y(4) 



(6.12e) 



Substituting Ap into ( |6.1 1| ) we have the momentum equation. For the D2Q9 model, 
as used by 



191, [TT|], h^^^ some simplifications, and finally 

.2 



dt{pUp) + daipUaU/s) 



dpp + vda [dp{puc) + docipup)] + F;3p + 



-7;9tipFfi) + [t - - ] dadt{pUaUp) - rdaipFfjUa + pFaUfs) (6.13) 



with 

c% ( 1 

u = r 

3 V 2 

The unwanted terms in ( |6.13[ ) are those of the second row, all of viscous order. 
Mach number M is the ratio of fluid velocity and speed of sound, which in Lattice 
Boltzmann models is proportional to lattice spacing c. In the method of successive 
approximations we have no expansion parameter e with which to scale Mach num- 
ber, but we can collect the terms in Aj3 with smallest Mach number, and therefore 
neglect the term in ( |6.13[ ) proportional to (r — ^). The momentum equations up to 
viscous terms hence read 



dt{pUp) + daipUaUp) 



-dpp + vda{pdpUa + pdaUp) + Fpp + 



I {-dtipFp) - rdo^ipFpu^ + pF^up)) + 0{5l) + 0{M^) (6.14) 



6.1 Summary 

• This model is Galilean consistent with small Mach number only. 

• Extra terms involving the force appear in the continuity and momentum equa- 
tions. To have correct Navier-Stokes equations we therefore have to assume 
that F is small, that is that we are close to the ideal gas. 
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7 An LBE scheme with high order accuracy 



The purpose of this section is to show that it suffices to slightly modify the Lattice 
Boltzmann model of section ^ to get correct viscous order equations. The model 
which we will describe was introduced in [ |T9| ] with a specific choice for the force 
term T . In [ ]TP| ] it is not stated very explicitly that these modifications solve the 
problem of unphysical viscous-order terms, and we therefore belabour this point 
here. 



We assume an isothermal situation. The starting point is to integrate equation ( |6.3[ ) 
to better than linear order for T 



/i(x + Qibt, t + 6t)- fi{^, t) -- 

•'t+St _ jeq 

A 



dt+ ^-^^^"V r^t 



(7.1) 



Applying the trapezoidal rule to the second integral of right-hand side (but not to 
the first), we have the Lattice Boltzmann model (compare (|6.4|)) 



5t 
+ - 
t 2 



RT 



fi 



eq 



F.(e^-u) 

t+5t RT 



(7.2) 



where r = A/(5t is the non-dimensional relaxation time. The right-hand side in- 
volves the quantities evaluated dXt + 5t. To eliminate this implicitness, we introduce 
the new variable: 



hi = fi 



F.{ei-u] 
2RT 



(7.3) 



in term of which the Lattice Boltzmann equation ( p.2[ ) is: 



/ii(x + ei6t, t + 6t) - hi{yi, t) 



hi^,t)-h-''i^,t) F.(e,-u) 
-r ^^J—^^ 



(7.4) 



Another interpretation of this procedure is that we evaluated all terms in equation 
( [7?2| ) to the same order, with a collision operator for the original variables fi that 
was actually not of the single relaxation time form. Instead it was something differ- 
ent that gave the single relaxation time form for the auxiliary variables h^. 
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The equilibrium distribution for h contains terms up to three velocities, and reads 

F.iei - u 



1-3- 



2c2 



-5t 



eq 



(7.5) 



where is the standard choice of (|3.6|). The macroscopic variables given by 
eqn.(r7T3D become 



m = ^ hiBi + ^pFSt 



(7.6) 
(7.7) 



We now wish to compute the macroscopic equations by Chapman-Enskog expan- 
sions. We will for notational ease use both hi and fi in the expansions. Since these 
two quantities differ by a term to order e (i.e. equations ( |7?3| ) and ([73|)), a choice 
has to be made for which of the two in which we count powers of e. The natural 
choice is to count in orders of e for fi, but to perform the expansion in hi, since the 
equations for hi are explicit. This will mix orders for hi, but in an unambiguous 
manner, as the expansion can be mapped back to order by order in e for /j. To 
order accuracy we have 

h, ^ hf^ + e/if ) + e'hf^ = hf-T [eVt + ^P?^ h, + er^. (7.8) 

where each term /i-"^ includes a term of order e, and where 

_ F.(e, -u) , 
RT 

Regrouping the right-hand side of eqn.(f7?^ by the leading order of e we have 

h^^ + e/^P + e'h^^ =hr - erV.ihf^ + ehf^ + e'hf^) 
-e^-^ihf^+eh^+^h^) 
+eTj^i (7.9) 



which is 



6°: 










= 






6^: 


= 


2 *o 


,,,0. 



(7.10a) 

rj^i (7.10b) 
rdt.hf^ - rVt.hf'^ (7.10c) 



37 



Equation ( [7.10a| ) here determines hf'^ to be of mixed zeroth and first order in e, i.e. 
given by ([7.5D. Since the hydrodynamic modes belong to to equilbrium distribution 
functions, we have for n greater than zero the constraints 



E =0 



(7.11) 
(7.12) 



For the lowest order hf^^ = h'^'^ we have on the other hand 



E''™^-=E^ 

i 



eq 

i ^ia 



lE^. 



--T^'^^p5ap + pUaU/3 - - 



pFocUp + pFpUa 



(7.13) 
(7.14) 



T(2) 



{F.u)pUaUp 



Summing equation ( |8.21| ), using equations d7.13D , (|7.14j ) and ( |6.5a| ), gives 



dtaP + daipUa) = -da{pFa 



(7.15) 



(7.16a) 



To second-order in e we get by summing ( p.lOc[ ) 



(0) 



(7.17a) 



The first and the third terms is substituted from equations ( |6.5aD , ( |6.5b| ) and ( |7.13D 
The second term is zero, and the two others give 



(7.17b) 



Combining the first and the second order results for hi by dt = dt^ + edt^ we have 
the continuity equation: 



dtp + V.{pu)=0 



(7.18) 
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valid up to order e^. In other words, there is no density diffusion term in this equa- 
tion. 

To first order, the momentum equation can be computed by multiplying equation 

(7.19a) 



2]\ i with ej/3, then summing over i, and using ( |7.15D : 



which is 



^^d^[{F.u)pu^up] + '-dt,{pFp) + pFp (7.19b) 



The second order momentum equation are obtained by multiplying equation ( p.lUcp 
with and summing over i, which leads to 



i i ^ ' i 

Substituting in equations ( p.l4[ ), ( |6.5b[ ) and ( |6.5cD , we have 



(7.20a) 



i(9t„(pF^) - (pFaUf^ + pFpUa - :^{F.u) pUaUp ) (7.20b) 



Combining the first and the second order results, equations ( [7.19bp and ( [7.20bD , for 

hi by dt = dto + edt^, we have, 

dtipup) + da{puaup) = -T(2)a^p + e - 1^ dji^^l + pFp (7.21) 
where n[^^j = CiaCiishf^ is the first-order momentum flux tensor. We have 

= E ^i«^«/3^i^^ ^ ~^ E ^ic^eiijVt^hf^ + ^ E ^io^em^i 

i i i 

= ~^ ( ^^0 E ^i«^«/3^f ^ + ^7 E ^i'y^iP^il^f^ ) + ^ E ^i^^iP^i 
\ i i / i 
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and in which the sums can be simplified to 

eiaeiphf^ = ^ eiaCifsff^ ~ eiaCipJ^i {1 .lid) 



^ ^ (^ia^^iP^i^hj ^ ^ (^ia^^ip^i^ f i 2 ^«a^*/3^i7*^i (7.22b) 



The second terms in equations ( |7.22p are order of e. Therefore, when we substitute 
them back into ( [7.21D , these terms become e^, and should be neglected. With this 



simplification n[^^j can be rewritten, 



n(^) = -r 



\ i i / i 



(0) 



where n^°j is zeroth-order momentum flux tensor, and the convective derivative has 
been computed above in ( |5.23p 

Aoni°j =dt,{puaup) + T(2) [d^ipu^) + d^ipua)] + 0{6t) (7.23) 
=T(2) [pd^u^ + pdpu^) + p{u^Fp + upF^) + 0{5t) + 0{v?) (7.24) 

The first-order momentum flux tensor is therefore in the small Mach number ap- 
proximation and to lowest order in e equal to 

= -rT(2) {pd^up + pdpuo) (7.25) 

Substituting back the above results into ( [7.21[ ) we have the momentum equation, 

dt{pup) + daipUaUfs) = -T^^^ dpp + vdaipdaUfi + pd^Uo) + pFp (1.16) 
where 

The final equation is the Navier-Stokes equation of an ideal gas with a force field. 
This force field can be an external force, but it can also be an internal force which 
describes capillary forces and non-ideal corrections. By comparison with (|2.12D we 



see that it suffices to choose pFp = T^'^^dpp + daTap to describe isothermal flow 
with reversible stress tensor Tap. If we for instance wish to have the capillary stress 
tensor of ( [2.10[ ) we should choose F = (T*^^) + (p(/ — p)' — KV'^p) V log p, where 



the Laplacian and the gradient can be computed by discretization on the grid as in 
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Summary 

This modifies second-order Lattice Boltzmann equation gives rise to correct 
continuity and momentum equations to viscous order. The local and nonlocal 
presure terms in a diffuse interface model can be introduced via the force field 
F. 

As in previous sections, a small Mach number expansion has been used, such 
that terms with three velocities are ignored. If desired this can be improved 
by using lattices with higher isotropy, following [^, [T^. 



The model is isothermal. 
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8 Thermal lattice Boltzmann 



8.1 General considerationsn 

The equations of non-isentropic fluid flow are the continuity equation, the momen- 
tum equation, an equation of state and either the energy or entropy transport equa- 
tions. In the context of diffuse interface theory, the energy and entropy transport 
equations are ([2.15|) and (|2.16D above. One difficulty in constructing a Lattice Boltz- 



mann scheme for non-isentropic flow is that we have so far not defined the internal 
energy or the entropy in terms of the distribution function. 

Let us start by remarking that temperature can be defined in terms of the thermal 
kinetic energy as 

P^T = fiiCia - uf (8.1) 

i 

where p is the mass density, and Et = ^RT. In the rest of this section we will as- 
sume units have been chosen such that the gas constant R is one. A first requirement 
of a nonisotropic Lattice Boltzmann model is that temperature is not a constant, i.e. 
that the equilibrium distribution function is modified to depend on T. 



In equilibrium, the internal energy and entropy densities per unit volume are re- 
lated by ([^, note that ej and s here stand for densities per unit mass) 

d{pei) = Td{ps) + pdp (8.2) 

The free energy density and the pressure on the other obey 

pf = p{ei-Ts) p=-p{f-p) (8.3) 

from which follows 

dei = -^dp + Tds df = -^dp- sdT (8.4) 
p^ p^ 

If we with slight abuse of notation denote the free energy functional density per unit 
mass in the Cahn-Hilliard theory by 

/[p,Vp,T] = /(p,T) + ^|Vp|^ (8.5) 
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where /(p, T) is the free energy density of a homogeneous phase at density p and 
temperature T, then 

P = V^|r = Vg-flVpr-iTpVV (8.6) 
^/[ ] , df 



6T dT 

We note that p is equal to the diagonal component of the stress tensor in the second 
formulation discussed in section |[ We also note that in equilibrium, as we assume 
here, the non-diagonal terms of the stress tensor are zero, and the stress tensor 
can therefore be considered a non-equilibrium generalization of pressure. Entropy 
density in diffuse interface theory is just a function of p and T, while internal energy 
density depends on density gradients, i.e. 

ej[p, Vp, T] = f[]+Ts = Slip, s) + ^|Vpp (8.8) 

where {ej, s) and (/, T) forms a Legendre transform pair: 

ej{p,s) = fip,T) + sT s = -^ (8.9) 

Thermodynamics in a system with mass motion follows from the free energy (see 
[§, section 8.4) 

F(r,V,M,v) = E-T^-P- V (8.10) 

where P is total momentum, and v is the velocity of the rest frame of the fluid 
relative to the laboratory. We have P = Mv, and the internal energy with total 
momentum P is 

E{S, V, M, P) = E{S, M, 0) + (8.1 1) 

from which follows 

F(T, V, M, v) = F(T, V, M, 0) - ^Mv^ (8.12) 

The internal energy ej, the free energy density / and momentum density per unit 
mass ^ in a system in motion are therefore related by 

/(p, T, v) = ei{p, s, 0-Ts- ei{p, s, = ei{p, s, 0) + (8.13) 



^In equilibrium the average momentum density per unit mass is the velocity v. 
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8.2 Non-isentropic Lattice Boltzmann 

The equilibrium distribution function can be generalized to depend on temperature 
TasFI 



1 + + fe,.u)2T(2) u2 



(8.14) 



T 2TTW 2T 

We consider the continuous-time Boltzmann equation in the form 

^ + i.Vf + ¥.VJ = nU) (8.15) 

where we identify ^ = d^e{p, s,^) and F = —dx£{p{x), s{x),^). The left-hand side 
of ( |8.15D is Liouville's equation. 

Here is a list of unanswered questions, that we feel is pertinent: 

• What quantities should the collision operator n conserve ? Presumably mass, 
momentum and energy ? 

• Can this then be done with a single relaxation time scheme ? 

• The internal energy density ei is a function of p and s, but s is a function of 
p and T. Hence we can consider si a function of p and T, both quantities 
which are naturally defined for the LBE. More explicitly, on the lattice, we 
could have 

K 

ei{x, t) = e%p{x, t),T{x, t)) + ^^^^ ^^^^^^ ^ Wi{p{x + e^, t) - p{x, t)f 

(8.16) 

where e']{p{x, t),T(x, t)) stands for the bulk contribution to the internal en- 
ergy, and the rest is a discrete approximation to the gradient terms. 

How to construct an LBE scheme that conserves such a beast? 



8.3 An attempt to generalize section 7 

The LBE models in the previous sections are all isothermal. The objective of the 
present section is to describe work towards adapting the model described in section 
to non-isothermal flow in a non-ideal fluid. The procedure follows the one for 

"^We assume units such that R = 1. 
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heat flow in an one-phase fluid with an ideal gas equation of state introduced by 
Alexander et al [jH]. This procedure has a limitation, in that the Prandtl number can 
not be varied. 



Several Lattice Boltzmann models of fluids with heat and mass flow have been 
introduced in the literature. Most of these share a defect of the one presented here, 
in that an equation is derived for the variation of temperature. The general equation 
of heat transfer ( |2.16D involves the variation of entropy density, and only reduces 
to an equation for temperature in the incompressible limit, see § 50. We will 
not attempt a review of these models, but refer to [ESI E^, EO, O, B^. The work 



of Palmer & Rector [ p4| ] deserves however to be singled out, as it uses a differ- 
ent approach where the Prandtl number can be changed, and simulates the energy 



transport equation ( |2.15| ), albeit without the effects of viscous heating and intersti- 
tial working. 



Let us start with defining the thermal kinetic energy as 



(8.17) 



where p is the mass density, and e = y-RT[|. The equilibrium distribution function 
only leads to isothermal flows, as e is then T^'^\ However, by changing the 
equilibrium function to be a function of p, u and T as in [jl|], we can have the 

^ As an aside, let us note that a natural definition of entropy density per unit mass would be 

i 

where in equilibrium s is a function of p and T. A Lattice Boltzmann scheme starting from such a 
definition would naturally lead to a convective derivative of entropy density, as on the left-hand side 
of ( 2.16| ). The right hand side of (2.16) however only follows if the equilibrium function is chosen 
in an suitable form: this work is not attempted here. 
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following set of constraints: 

E/r=P (8.18a) 

i 

J2f!'eia=puc, (8.18b) 

i 

X] fi'^^ic.eip =pUcUp + pRT5af3 (8.18c) 

i 

f^eiaeipei^ =pUaUpU^ + pRT{6al3U^ + Sa-yUi3 + 6i3.yUa) (8.18d) 

i 

Ereq 2 2 , 

i 

(D + 2)p{RTf5p^ + p(i?r)M^5^^ + (D + 4)p(i?T)M^M, 

(8.18e) 

We now introduce an auxiliary quantity hi as in section ^ which obeys 

Y.hT=p (8.19a) 

i 

K'e,^ =pu^ - ^F^p{RT)6t (8.19b) 

i 

^ h1''eic,eip =pUo,up + pRTS^/s - -^^{F/sUa + Fo,up)pRT5t (8.19c) 

and a force term in the Lattice Boltzmann equations ( |7.4| ) which obeys 

^J^i=0 (8.20a) 

j 

5^-^^ei^=;^pi?rF^ (8.20b) 

i 

^ TitidCi^ =^^pRT{FfiU^ + F^M^) (8.20c) 
5^ T^e\^^^,i {{D + 2)p{RTfFp + p/^Tw^F^ + 2pRTF^u^up)) (8.20d) 

i 

We compute the macroscopic equations with Chapman-Enskog expansions as in 
section ^ The constraints we use to close the equations are the zeroth, first and 
second moments of the collision operator, that is 

• - f^') = - /^')^- = - /r)e..e,„ = (8.21) 
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The first and second order in expansion parameter e of the continuity equation read: 



dt,p + d^{pu^) = t^^^d^{F^pRT) (8.22a) 
dt,P = -^d^iF^pRT) (8.22b) 

Combining the first and the second order results by dt = dt^ + <^dt^ we have the 
continuity equation: 

dtp + V.ipu)=Q (8.22c) 
The first-order in e of the momentum equation is: 

dtoipup) + da{pUaUp) = - dpPo + ^i^^{FpUa + FaUp)pRT + 

^^dt,{FppRT) + ^FppRT (8.23) 

with the ideal part of the pressure Pq = pRT. The first-order in e of the temperature 
equation is calculated by taking the second moment of ( |8.21| ). 

•9*0 (^pRT + ^pul 



- (^pulup + ^^-^^RTup^ + ^pRTFpUp + J A„ ^ Ticl 



2 

ia 



(8.24a) 

To second-order, the temperature equation can be computed by the same way from 
equation ( |7.10cD . 

dt, f^pRT + ^pul 



1 / 1 
2 



- E ^S'^^- - \l^to E ^^^l + E -^^^^ (8.25a) 

^ ^ i i i 

Combining the first and the second order of energy equations ( |8.24aD ,( |8.25aD by 

dt = dto +edt„we have, 

dt {^pRT + ]^pul^ = - ^dp (^pulup + ^^^±^RTu^ + 
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In the equation above, a term with the second-order of e: ^dt-^ J2i ^i^la 
glected. By using the first order of continuity and momentum equations ( |8.22aD , 
( |8.23[ ) we have 



(8.27) 

where P stands for the full effective pressure given by daP = da (pRT) — ^^27 pRTFa . 
Subtituting equation (|8.27p back into (|8.26|), we have, Pi 



~2 



dti-pRT) + dp{peup) = -Pdpup + 1 ("1 _ 1 J 2^ Ki'eia + 



^peFpUp (8.28) 



With the assumtion that, the local energy is conservative, the second term in r.h.s of 
equation above can be rewriten as, 

i \ i i i / 



The high-order terms of e can be neglected as same as ( [7. 22] ), 



\ 



\ 



dll 



(3) 



I3~t 



Some terms in the equation above can be further simplified as: 

dtipulup) =UaUpdt{pUa) + UaUp[dt{pUa) - UaOtp] + - Ufsdtp] 

= - d.,{pului3U^) - 2UaUi3daP - U^df^P + 



12 ^ 6 Or. 

-jj-peuaU/sFa + ^peu^Fg 



(8.29) 



From here on the calculations have not been rechecked thoroughly. The reader beware that the 
likelihood of accidental is not negligable, and, in particular, that there may be some confusion below 
in the use of RT and e. 
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dt{peu/3) =updt{pe) + e[dt{pup) - updtp] 

= - da{peUaUp) - ^peUpdaUa 



D 



edfsipe) + 



-ipprpeFaUaUfS + -ipfzpS FfS 



(8.30) 



By using the equations (8.29), (8.30) and (8.20d), we have, 



i ^ia 



--^peupdaUc 



4(D + 2) 
6 



ped^e 



GjD + 4) 2 



(8.31) 



Subtituting back equation ( |8.31D into ( |8.28| ), we have the enery equation, [| 

dt{pe) + d/sipeup) = -Pdpup + dp{Kdpe) + dp[p{di}Ua + daUp)ua\ + 



dfi{XupdaUa) + -(/i - \)df3{eFp) + -pdp{FaUaUp) + —peFpUp (8.33) 



where. 



4 

2f2 + D 



(8.34) 
(8.35) 
(8.36) 



^ Footnote added: the equation of heat transfer in an incompressible fluid is (compare [p4(], eq. 
50.2) 



pDtT = {k/cp)AT + (/^/2cp) {daV/s + dpVaY 



(8.32) 



where Cp is the specific heat at constant pressure. Equation (8.31) definitely contains more terms 



than this. Among these are terms of the right hand side of the energy equation (2.15), but also 
other terms. The meaning and validity of equation ( |8.31 ) is unclear. The calculations of this section 
however well illustrate that constracting thermal multiphase LBE models is not trivial. 



49 



9 Optimizing a diffuse interface model 



We have seen that a diffuse interface theory of a multiphase flow involves a free 
energy density, which in the isothermal situation is a function of density. For a real 
fluid this energy density is an observable quantitity in a homogeneous phase. In the 
interface, the real form of the free energy density is less constrained by experimen- 
tal data, and one has some freedom of choice in the model. 

The numerical difficulty in solving a problem of a multiphase flow with an interface 
by a Lattice Boltzmann, or other diffuse interface method, depends on the interface 
thickness ^, since the grid size Ax cannot be less than ^ 0. A typical density profile 
in an interface is shown in Fig. [T|. The interface thickness is approximatively related 




' ^ * normal direction (z) 

Figure 1 : Density profile 

to the density gradient and the density difference between the two phases as 

dp Ph - Pi 
dz ^ 

In a typical application {ph — pi) is something we would like to keep fixed. Simi- 
larly, surface tension, and the physical properties of the two phases, are also things 
we would like to keep fixed. 



(9.1) 



From ( p^ and ( [2.18D (or ( |2.22D ) follow that the surface tension can be written 

2 j\pf-pii + P)dz (9.2) 



a 



**We assume a homogeneous grid, although that is not necessary in the Lattice Boltzmann method. 
An adaptive grid would involve some sort of front-tracking, and lead to problems outside the scope 



of this paper. 
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where 2 is a coordinate in the normal direction to an interface. The thermodynamic 
rule that the pressure and the chemical potential in two coexisting phases are equal, 
implies that the two phases (at densities pi and p2 and bulk free energy densities ipi 
and 11)2) are related by the the double tangent construction of Fig.^. Let us introduce 




Density 



Figure 2: A typical bulk free energy for a one component system 

-0, a linear interpolation of if), and write 

^(p)=^(p) + A^ 7/i = ^i+p(p-pi) (9.3) 

With P the pressure in one of the homogeneous phases (for instance the 1 phase) 
we have in the interface region 



= ip - Pfi + P 



(9.4) 



Equations ( P.4D and ( |9.2D express that the excess free energy per unity of surface 
area of the interface is the excess free energy density in the interface region, inte- 
grated over the normal direction to the surface. 



By using ( p7T| ) we have 



(9.5) 



where A-^ is the mean excess free energy density in the interface region. Equation 
(Oh on the other hand means 



Alp K 



(9.6) 
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If the parameters to be kept fixed are a and Ap, the one to be taken large is ^, and 
the two parameters to be freely adjusted are K and A^/^ we have 



a 



Alp - K 



(Apy 



(9.7) 



To have ^ large, all else constant, one should therefore take Atp small and K large 
in the same proportion. The optimization problem addressed here is that of maxi- 
mizing the density difference Ap, while keeping a and ^ constant. Aip should then 
be held constant, while K should be inversely proportional to the square of Ap. 

If we use a standard equation of state, such as the Van der Waals 

P 



P 



^-ap" ^ = pRT\n 

1 — 



1 - bp 



ap 



(9.8) 



the interface is thin, except close to the critical point, where however also the density 
difference between the two phases is small. Let us now consider the following 
"stretched" Van der Waals free energy, where the free energy curve is a straight line 
in the interval [pm — Ap, pm + Ap], see Fig.^ 




Figure 3: A modification of bulk free energy for a one component system 



Ipnewip) = 1p{P + Ap) p<Pm-Ap 

Pm ) Pm 

Ap<p<pm + Ap (9.9) 

i^newip) = i'ip - Ap) p>Pm + Ap 

Here pm is the point of maximum excess free energy in the Van der Waals energy, 
and A'ip(pm) is that excess energy. The effect of the perturbation is hence to move 
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apart the two densities of the coexisting phases. 



The pressure in the stretched region is constant, and equal to the pressure in the 
two coexisting phases, see Fig.H With this method, the density difference betweem 



-d' 




Density 



Figure 4: A modification of equation of state 
two phases can be changed quite easily, see FigJ| 



Origin 




Figure 5: Solutions with the different modification of Ap = 0, 1, 2 



9.1 Summary 

The form of the free energy in the interface determines surface tension and interface 
thickness. In a physical liquid these two quantities vary together with temperature 
and pressure in the phase diagram, and the interface is only thick close to the criti- 
cal point. In a numerical procedure that demands the resolution of the interface, the 
interface thickness sets a limit on the smallest grid size that can be used. 
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The free energy in the interface can be modified, so as to e.g. keep surface ten- 
sion and interface thickness constant, while increasing the density difference. From 
a numerical point of view, this procedure amounts to smearing the interface, keep- 
ing all else constant, i.e., in some sense the opposite approach to the use of singular 
interface theory. The advantages of the procedure are that problems which are dif- 
ficult in the singular interface theory can be readily simulated (see list on page 4), 
in particular there is no difficulty with changes in topology. The drawback is that 
the simulations will be of a fictitious liquid with thick interfaces, and the validity 
of the simulation results would in the end have to be established from comparison 
with experiments. 
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10 List of symbols and conventions used 





speed of sound 


Gi 


lattice (particle) discrete velocities 


fi 


discrete one-particle distribution function 


P P 


pressure and pressure tensor 


R 


gas constant 


T 


temperature 


u 


macroscopic velocity 


P 


density 


£i 


internal energy 


£k 


kinetic energy 


e 


expansion parameter 


0{...) 


on the order of ... 


Wi 


weight of different sub-lattices 


a, (5 


space coordinates 




Kroenecker delta symbols 


5t 


time scale 


Sx 


length scale 


V 


first viscosity 


K 


elastic stress coefficient 


X 


bulk viscosity 


V 


kinetic viscosity 


i 


second viscosity 


a 


surface tension 


T 


relaxation time 




bulk free-energy density 




non-local free-energy functional 




collision operator 


F 


force term 


n 


momentum flux tensor 


T(2) 


Lattice constant, second order tensors 


Y(4) 


Lattice constant, fourth order tensors 



Greek indices (a,/?, . . .) generally label spatial coordinates, while Latin indices 
(i, j,. . .) label lattice vectors. The Einstein convention of summing repeated spatial 
indices has been used throughout the paper. Hence, if a. — {ai, . . . ^an} and b = 
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{bi, . . . , bo} are two vectors, then Uaba = a • b = Z]q=i ciaba- Repeated indices 
labeling lattice coordinates are not summed over, unless so indicated. 
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